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ABSTRACT 


This report describes the research performed during the period March 
1987 - July 1988 under a joint research program between the Mechanical 
Engineering Department of the University of Maryland and the Center for Fire 
Research of the National Bureau of Standards. The research is conducted in 
the laboratories of the CFR by a Graduate Research Assistant of the ME 
Department under the joint supervision of Dr. Marino di Marzo (ME Dept. - 


UMCP) and of Dr. David D. Evans (CFR - NBS). 


The formulation of a model for the prediction of the cooling induced by 
an evaporating droplet impinging a semi-infinite solid is the subject of this 
report. The thermal interactions during the evaporation of a liquid droplet 
deposited on a low conductivity semi-infinite solid are complex because the 
evaporative process is coupled to the solid intense local cooling. Numerical 
techniques based on finite difference methods have failed to provide 
meaningful results. This is due to the sharp temperature gradients in the 
proximity of the droplet edge which cause instabilities in the solution for 
reasonable time steps due to the explicit coupling of the liquid-vapor 
regions. An integral method was proposed by Dr. Baum (CFR - NBS) in order to 
overcome these difficulties. The methodology and its application to this 
specific problem is described in detail. Preliminary result will also 


illustrate its validity. 


A brief note on the convective heat transfer coefficient measured in the 


previous experiments on aluminum and Macor is also included in this report. 
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iy INTRODUCTION 


aire’ Background 


The long term objective of the study of droplet-solid interaction is to 
obtain information applicable to the extinguishment of fire through a droplet 
array (e.g. spray). The solids of concern include porous materials, typical 


of fire applications. 


This report describes the research performed during the period March 
1987 - July 1988 under the joint research program conducted by the Center for 
Fire Research (NBS) and the Mechanical Engineering Department (University of 
Maryland) addressing the Transient Cooling of a Hot Surface by Droplet 


Evaporation. This joint research program was initiated in January 1985. 


The research described hereafter constitutes the first phase of an 
extensive research program aimed at developing accurate droplet cooling models 
of burning solid fuel surfaces. Many studies have been performed to quantify 
the vaporization process of both single and multi droplet arrays impacting on 
hot surfaces. For the studies found in the published literature, the full 
span of the droplet vaporization processes are usually reported. These would 
include evaporation, nucleate boiling, film boiling and te laentrose 
transition. This present research is more limited in the span of vaporization 
process studied, being only concerned with the evaporation of single droplets 


on a hot surface. 


» 


Several important results were obtained in the first two years of 
research. In particular, the modelling of the boundary condition at the 
liquid-vapor interface (at the droplet exposed surface) was validated with the 
data collected for water droplets evaporating on an aluminum block (diMarzo 
1986a, 1986b, 1988). The simple model that describes the cooling effect 
induced in the aluminum by the evaporating droplet (diMarzo 1986a, 1987b, 
1987c) provided the input for the formulation of a preliminary multi-droplet 
model (diMarzo 1987a). An extensive experimental data base was also collected 
for water droplet evaporating on Macor (a glass like material able to 
withstand strong thermal stresses). These data will be used to validate the 
more complex model required to describe the coupled liquid-solid interaction 


for the Macor case. 


1.2 Qverview 

The rational for the more complex coupled model has been described at 
length in the previous reports. The low thermal conductivity of Macor 
invalidates the assumption of constant uniform temperature under the liquid 
droplet which was crucial in developing the simple model for the case of 
aluminum. This fact requires that the droplet evaporative process be solved 


jointly with the cooling transient in the solid. 


Two unsuccessful attempts to use finite difference schemes for the 
solution of the coupled droplet-solid problem were performed during the later 
part of 1987. The first numerical technique, based on an Alternate Direction 
Iteration (ADI) scheme, was implicit both in the liquid and in the vapor 


region, but required an explicit coupling at the liquid-solid interface. The 


second numerical technique was also implicit and based on the iterative Gauss 
Siedel (GS) scheme applied to the coupled liquid and solid field. Both 
techniques failed because the sharp localized gradients at the droplet edge 
were introducing unmanageable oscillations for reasonable time steps ( = one 


second). 


A different approach based on integral methods was proposed by Dr. Baum 
to overcome these difficulties. Two important advantages can be identified in 
this methodology: a) only the surface thermal behavior are considered in the 
solution; b) the integral nature of the technique is better suited to handle 
sharp gradients which caused oscillatory behavior when differentiated. This 
integral method is based on the introduction of appropriate Green's function 
to transform the problem from spatial to superficial via the Gauss theorem. 
Once the problem is cast in superficial terms, then it reduces to a double 
integral in time and over a closed surface of an argument composed of the 


original variable and of the Green’s function. 


A brief note on the convective heat transfer coefficient measured in the 
aluminum and Macor experiments opens the report. This note is necessary for 
completion and to summarize the various results concerning this measurements 
which is fundamental for-the derivation of the mass transfer coefficient at 


the liquid-vapor interface (at the droplet exposed surface). 


1.3 Program Development 
In this section a brief synopsis of the various research activities is 


given in the time frame of the reporting period. 


Period March 1987 - August 1987 Derivation and implementation of the ADI 


scheme with repeated unsuccessful attempts to obtain converging solutions with 


acceptable time steps. 


Period August 1987 - December 1987 Derivation and implementation of the GS 


scheme leading also to unsuccessful results. 


Period January 1988 - July 1988 Derivation and implementation of a boundary 


element method via an appropriately defined Green’s function which is 
successfully validated by existing data obtained for the aluminum experiment. 
Development of a short note comparing the experimental data for the convective 


heat transfer coefficient. 


Ge OVERALL AND CONVECTIVE HEAT TRANSFER COEFFICIENT 


A Overall Heat Transfer Coefficient 


The overall heat transfer coefficient is evaluated by measuring the heat 
flux trough the semi-infinite solid and the solid-surface-to-ambient-air 
temperature difference. The details of these experimental measurements are 


described in previous reports (diMarzo et al. 1986a and 1987a). 


The objective of this brief summary of these results is to validate the 
data concerning the convective heat transfer coefficient measured in two 
different experimental apparata (aluminum and Macor) by different individuals 
in different time frames. This validation is particularly important because 
the theoretical predictions available in the literature (McAdams 1957) are 


consistently underpredicting these data. 


In order to deduce the convective heat transfer coefficient from the 
overall heat transfer coefficient data it is important to recognize that the 
heat is removed from the solid surface by convection and radiation in 


accordance to the following relationship: 
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where the emissivity of the surface ‘ce’ is 0.17 and 0.84 for aluminum and 
Macor respectively. This values were obtained using the Gier Dunkle Infrared 


Reflectometer. 


Zz Experimental Measurements 


Equation (1) allows to calculate the convective and radiative heat 
transfer coefficients for aluminum and Macor. The radiative heat transfer 


coefficient is defined as: 


hee oe wencreeT eT. + T.) (2) 


r 
The results are presented in Tables 1 and 2 respectively. 


Lee Theory and Data Comparison 


Theoretical values of the convective heat transfer coefficient for the 
aluminum block were found by considering the exposed surface as a disk. This 
can be done since the sides of the block a ee insulated and the heat loss 
from there is negligible. A correlation for the convective heat transfer 
coefficient in natural convection from a hot disk facing upward is proposed by 


McAdams 1954: 


T. 0.25 


ye ghey (—3P ? (3) 


where S/P is the ratio of the disk surface and perimeter. The comparison of 
the data for Macor and aluminum with this correlation is presented in Fig. 1. 
The theoretical results are consistently lower than the experimental results 


probably due to slight laboratory air circulation. 


Table 1 - Aluminum Heat Transfer Coefficients 
To) 25°C e= 0.17 
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Table 2 - Macor Heat Transfer Coefficients 
Tan=323°C e= 0.84 


CONVECTIVE HEAT TRANSFER COEFFICIENT 


THEORETICAL VS EXPERIMENTAL 
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FIGURE ONE - Convective Heat Transfer Coefficient: Theoretical Versus 
Experimental (— theory, ® aluminum data, & Macor data) 
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a MODEL OF DROPLET SOLID INTERACTIONS 


a The Governing Equations 


The conduction equation describes the thermal transient in a solid and 
in a motionless liquid. The full description of the problem requires the 
identification of all the boundary conditions and the coupling of the liquid 
and solid regions at the wetted surface under the droplet (see Fig. 2). The 
conduction equation can be written for both the liquid and the solid, keeping 


in mind that the thermal diffusivity should be changed accordingly, that is: 


One 2 
anew oC ov OL (4) 


The liquid region is bounded by a liquid-vapor interface where the mass 
transfer is coupled with the heat transfer as extensively described in the 


previous reports. This boundary condition can be summarized as: 


Dee he, Xie 
pee a (eo ee 0 2624) he Stor co, ina (5) 


The liquid-solid boundary conditions express the continuity of the 


temperature and the conservation of energy across the interface, namely: 


is Ne ae (6) 


TT (7) 


FIGURE TWO - The geometry and its coordinates 
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The boundary condition at the solid-air interface accounts for the 


convective and radiative heat transfer contributions: 


SAK SENT SERRA TAS, Teta er(T* ef Tz) (8) 


The axis-symmetric nature of the problem grants that, both in spherical 
and cylindrical coordinates, the gradient of T is zero on the vertical axis 


through the origin of the coordinate systen. 


The initial condition can be a linear one-dimensional temperature 
distribution in the solid and uniform temperature in the liquid or uniform 


temperatures in both liquid and solid. 


The droplet induced cooling of a semi-infinite body has a dual time- 
spatial parametric description. One could select as characteristic length 
some length related to the droplet wetting surface or the ratio of the solid 
thermal conductivity to the convective heat transfer coefficient solid-to-air 
or even some mixed defined length such as the so called penetration length (1 
= [a t]°->) where the time could be the droplet evaporation time and a is the 
solid thermal diffusivity. Similarly the time scale could be the evaporation 


time or a grouping of solid related properties such as [(p c, k)/h?] or again 


P 
a mixed definition such as the square of some length related to the droplet 
wetting surface over the solid thermal diffusivity. The point of all this 


discussion been that non-dimensionalization of the governing equation is not 


univocal and it should not be pursued at this preliminary stage. 


ll 


3.2 The Solution Scheme 

The difficulties encountered in the application of differential methods 
to the solution of the coupled solid-liquid thermal transient suggest that a 
different approach should be considered. Two major considerations lead to the 
selection of the appropriate methodology: a) the relevant event are taking 
place on the surface of the domain; b) the nature of the problem with its 


sharp localized changes in the thermal gradient requires an integral approach. 


The idea of reducing the problem to a surface problem is appealing 
because only the points on the solid (and liquid) surface can be monitored 
experimentally. An infrared thermographic technique will be used to monitor 
the surface temperature during the droplet evaporation transient. This 
technique is non-intrusive and is highly suitable to this particular problem. 
Therefore, the computation limited to the surface points is more efficient and 
allows a more precise definition of the noding in the region of concern (eg. 


at the droplet outer edge). 


An integral methodology is also desirable because the solution at a 
given point (eg. point temperature) is obtained by superimposing a large 
number of contributions from the neighboring points hence the localized, 
drastic thermal changes of this cooling process are smoothed out. When a 
finite difference technique fe used the sharp gradients are locally amplified 


thus causing the observed instabilities in the solution. 


To obtain the desired result, it is necessary to introduce an adjoint 


equation in terms of the Green’s function ’'G’, that is: 
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5 =-avV G (9) 


By multiplying Eq. (4) by G and Eq. (9) by T respectively and integrating over 


the domain making use of Gauss’ theorem, one obtains: 


$02 a aca] five - oo + oes a (10) 
ie. ANY t js 


At this point the left-hand-side is expanded into two volume integrals at time 
t and at time zero respectively. By G at time t and on the boundary as the 
Dirac function, the volume integral at time t reduces to T. Further, one can 
introduce a new variable u in at of the temperature which is zero in all the 


volume at the initial time (t = 0). For instance one can define u as: 
Meee gq 2/7 kK) Gil) 


With these two modifications the final result is achieved in the form: 
ua te *. (G Vu - u VG) e nds dt C12) 


The solution of the transient conduction equation in a spatial domain of 
volume ‘V’' has now been reduced to an integration over the surface ‘'S’ that 
bounds the domain. The selection of an appropriate function 'G’ will be 


ps: 


discussed in the following. However, consider that the function ‘G’ must 
approach the Dirac function when the time is the present time. For the time 


being let us say that one possible selection is: 


G(¥,t:¥,,.t,)-(4"at,) be | 4a t, (13) 


where the vector position ‘Vy’ is in terms of any specified coordinate system. 
Note that G depends on the point v which is the point of concern and on the 
point v, which is the point identified by the integration dummy coordinates. 
Further, in order for Eq. (9) to be satisfied a new time scale t, must be 
introduced. The origin of this time scale is at the present time t and t, 
increases in the negative directions of time. Due to this set-up, t, is 
identified as the recollection time. Its significance will become apparent in 


the following. The dimensions of G is [meters]~* which is consistent with the 


dimensions required by Eq. (12). 


os] Boundary Conditions 


The surfaces of concern in the problem are identified in Fig. 3, namely 
they will be referred to as the cap (liquid-vapor interface), the disk 


(liquid-solid interface) and the plane (solid-vapor interface). 


The boundary conditions at exposed surfaces will be introduced in the 
formulation of the terms u and Vu in Eq. (12). The liquid-solid coupling 
conditions will be used implicitly in the solution scheme. This procedure 


will be explained in the following. 
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FIGURE THREE - 


Domains and surfaces 
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3.3.1 Solid-vapor interface This surface is a plane and a disk of radius ‘R’ 
has been cut out from it under the droplet. The appropriate coordinate system 
to integrate Eq. (12) is cylindrical, as dicteted by the problem’s geometry. 


A significant simplification can be introduced by selecting the G function as: 


2 
[ (xX - Kear (y - y,)°+ (Z. =< 'ZR8 ] 
G(¥ trv, C.J = (4 mane neat on 4a t, 


es - x)" + (y - ee (z + =)" 
+e. 


4at, (14) 


The clear advantage of this selection is that at z = 0 the term VG is zero. 


This means that Eq. (12) reduces to: 
u=a Ir <6 Vu e nds dt GP Bey 


This integral is not performed over a closed surface. The effect of points on 
far away from the droplet is negligible. Therefore, the surface bounding the 
domain in the far field is not included or its contribution is considered nil. 
{iio See the integral is performed over the surface 's’ which excludes the 
disk under the droplet. This portion of the solid surface will be considered 
with its adjacent liquid surface in the next section. With the implementation 
of the cylindrical coordinate system the equation is further simplified since 
the angular integration is reduced to a Bessel function term. The final form 


results in a double integral in time and along the radius ‘'r’ which can be 
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casted in the following form: 


t 
u(r,t) = os [ vuce, ,»t,) Ty pabhue 
J/4na eR 
2 
(a2 es cine) ] 
te] ee 
L, pesreate é 4at, dr, dto (16) 


There are singularities at t, = 0 and for r=r,. These singularities will be 


dealt with in the following. 


3.2.2 Solid-liquid interface The solid-liquid interface is under the base of 
the droplet. Formally Eq. (16) fully describes these points but the limits of 


integration should be modified as: 


t pR 
foie = roe 
J/4na 0°0 
2 
2x rx, gees 48) 
L, ATTEN e 4 at, dr, dto Glis) 


However, it is clear that both the temperature and its gradient are not 
known. By introducing the boundary conditions, Eqs. (6) and (7), we reduce 
the number of unknowns and we can solve the complete set of equations for 
temperatures and temperature gradients. The details of this procedure will 


become apparent in the following. 


For the solid the surface that bounds the domain is closed far away from 


jy 


the droplet. This means that the effect of points located on such a closing 
surface is negligible. The points on the disk and on the plane are the only 
one able to contribute to the solution. This is not the case for the liquid 
where the surface that bounds the domain is finite and where all the points on 
the this closed surface contribute to the solution. The impact of points 
belonging to one interface on points located on another interface will be 


discussed later on. 


Of particular interest is the nodalization at the corner of the droplet 
where liquid, solid and vapor meet. This zone will be handled as described in 
Fig. 4. Therefore each node accounts only for the portion of the surface to 


which it belongs. 


3.2.3 Liquid-vapor interface This surface is a portion of a sphere and both 
G and its gradient are non zero on this surface. Therefore it is appropriate 


to use spherical coordinates, which applied to Eq. (12), yield: 


2 . 2 os [cos(¢? - $,) - 1] 
u(¢,t) = =| {rue & 4a t, 
oy wae 0 40 
2 
2 Pe [COS(H™- oO.) eal) 2 ae (cos¢ cos¢, - 1) 
ve ETT lecnacee ow) hat . 
; ° 
2p, [cos(¢ - 4) - 1] 2 Be (cos¢ cos¢, + 1) 
ont amt hans, Sida Ws thu sp tect ‘ea 
2 
2 p, [cos(¢? - ¢,) +1] moe 2 po sing sing, 
e Gat, t oA 4at datings Gaadts 
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FIGURE FOUR - Droplet Edge Nodalization Detail 
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This formulation is far more complex than the previous ones due to the double 
term in the original equation. This equation is implicit in the unknown ‘u’ 
and the boundary condition is a functional relationship that calculates Vu as 
a function of u as it can be seen from Eq. (15). In the following, a 
conceptual simplification of all these relationships is sought in order to 


derive a simple algorithm to solve the problem. 


3.4 Weight and Forcing Function 


In order to simplify the task of handling these complex surface 
integrals, a decomposition in two portions is proposed: the value of the 
forcing function (be it u or Vu) is assumed to be constant in the small 
intervals dr,, dé@, and dt,. This approximation is valid in the spirit of 
finite difference approximation of the integration process. The various 


surface integrals can be recast in the form: 
n 
bie yueF aw oF, (19) 
i 


= 
where ‘W’ is a weight matrix and 'f’ is the vector of the forcing functions. 
It should be noted that for the liquid-vapor interface the formulation is more 


complex because double terms will be present due to the double forcing 


functions, namely u and Vu. The summation term is known since it involves 
previously calculated parameters. The second term on the right-hand-side 
contains unknowns. Note that the values at the liquid-solid interface are 


unknown as well as all the forcing functions on external surfaces which are 
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dependent on the value of the surface temperature itself. However the 
important simplification is in that once the weight are assigned at each pair 


of points they do not change at each time throughout the computation. 


In light of the previous observations it is important at this point to 
gain some insight in the concept of recollection time. The real time is 
elapsing as the evaporative cooling process takes place. The recollection 
time is zero at the present time and stretches its positive axis in the past. 
Further, the effect of past forcing functions fades as real time goes by, 
hence the weights of the past forcing function keep decreasing. Depending on 
the type of material and on the extent of the cooling, the recollection 
(memory) of the system is longer or shorter but in any case limited only to a 
finite number of previous time frames. Figure 5 clarifies the two different 


time coordinates. 


In general, each solution at a given point depends on the effects of all 
the other surface points. So far the effect of points belonging to the same 
surface have been considered. The inter-surface dependance will require 
special treatment as described later. More immediate consideration must be 
given to peculiar situations: The effect of the forcing function at the point 
of concern and the singular behavior of points at the origin of the spacial 


and temporal coordinates. 


3.4.1 r= xr _and ¢ = ¢ points The first particular case of interest is when 


the point of concern coincides with the source point. This corresponds to the 


maximum influence of a source point because the distance is minimal. In order 
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to treat this singularity in detail, an analytical integration of the general 
form of 'G’ and ‘VG’ is performed over the elementary surface portion 
represented by one node (eg. x - Ax/2 <= x = x + Ax/2). For the disk and the 


plane the result is: 


W(rer,,t,) = r ceadran (20) 


Note that the term At, is consistent with the formulation of Eq. (19) where 
the forcing function f is multiplied by the weight W. By comparing Eq. (19) 
with the general formulations for u, it is clear that a temporal and spacial 
integration must be performed on the weighing function. Such integration in 
this case is performed analytically in the space and numerically in time. The 
ne of the weight W are those of length as expected. In a analogous way 
the weight functions for the spherical cap are obtained. For this case the 


fact that ¢ = ¢, yields some readily noticeable simplifications in Eq. (18). 


3.4.2 t= 0 points and axis of symmetry points These singularities are 


addressed by staggering the discretization of the domain. The region 
neighboring the axis of symmetry is handled by locating the node nearest the 


symmetry axis is located at Ar,/2 or A¢,/2. 


Similarly the time discretization is not carried to the time origin (eg. 
the present time or t, = 0). In the following sections a detailed description 
of the temporal behavior of the weight function will be presented. In 
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particular analytical and numerical integrations over the initial time step 


will be compared in light of specific approximations of general formulation. 


3.4.3 The inter-surfaces dependance Going back to Eq. (12), the integral on 


the right hand side must be over a closed surface. This means that the effect 
of points located on the cap must be accounted for when computing the 
temperature of pOlnesten the liquid disk. 

In order to evaluate these contributions it is necessary to express the 
weight functions when the point of concern does not belong to the surface 
inherent to the coordinate system. As an example, consider the weight 
associated with a point on the liquid disk as influenced by a point source on 
the cap. In this instance the coordinate system must be spherical (since the 
optree point is on the cap) and the point of concern must be identified in 
spherical coordinates too. The poreace where the point source is located 
determines the coordinate system and all the previous formulation must be 
generalized to account for these points. In particular, for a point on the 


spherical cap and a source on the liquid disk the weight function is given by: 


Zz 
Zz ‘ei a abe 
65 54 ante 1.5 Corer $ ° 
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The distance from the disk to the point on the cap is z and it shows in the 


first term of the right hand side which was not present in the previous Eq. 
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(17) where the point of concern belongs to the disk (z = 0). 


3.5 Solution of the Equations 


In order to solve the liquid-solid coupled problem, the interfacial 
points must be related through the appropriate boundary conditions. This 
means that while the function u is explicitly known the plane, its values on 
the two coincident disks is linked to the value of Vu which is also unknown. 
The implicit nature of the solution on the cap is already apparent from 


inspection of Eq. (18). 


By rearranging the equations for all the points, ome can recast the 


problem in the form 


(22) 


where A is the modified weight matrix at the present time. This matrix is 
composed of all the weights associated with the unknown fluxes and with 


portion of the identity matrix associated with the unknown temperatures. 


The vector B is a vector that encompasses known quantities at present 
and past times as shown in the summation term of Eq. (19). Further B also 
includes the known weight-flux products or the known temperatures at the 
exposed surfaces (cap and plane) which are described by the boundary 
conditions. An example of this solution scheme is given in the following 


section. 
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4. PRELIMINARY RESULTS 


4.1 Weight Function Characterization 


The formulation of the problem in terms of weighing functions and 
forcing functions is expressed in Eq. (19). It is now necessary to carefully 
consider the behavior of the weighing function in order to itemize its spatial 
and temporal characteristics. ‘Ie is intuitive that the influence of a source 
on the point of concern is proportional to the distance and also inversely 
proportional to the elapsed recollection time. As the recollection time 
becomes larger the effect of more remote sources is felt but the influence of 
each individual source becomes less important. In the limit when the 
recollection time goes to zero (the now time), only the source at the point of 


concern is felt. 


4.1.1 Composition The weight function for the plane and the lLiquid-solid 
disks is composed of four terms: 
a) The recollection time t, to the -3/2 power. The exponent 3/2 is 
due to the chosen three dimensional Guassian form of the weight. 
b) The radial distance. 
c) The Bessel function L,. This function is the product of I, and 
the exponential of the negative argument. 
d) The exponential related to the Gaussian form where distance and 
time are ratioed. This term accounts for the distance between the 
source point and point of concern. 


The Bessel function c) and the radial distance b) account for the volume 


encompassed in the cylindrical coordinates. 


4.1.2 Temporal behavior The time characteristics which is of primary interest 
is presented for the most elementary combination of the above time dependent 
terms in Fig. 6. It is apparent from the graphical presentation that the most 
critical behavior (significant contribution) of the weight function occurs at 
the initiation which basically translates into the first time step. Therefore, 
the weight function has to be meticulously evaluated in the neighborhood of t, 
= 0. It is important to note that the weight function will converge to zero as 
the recollection time tends to zero. This means that the weight function is 


bounded throughout, which is a required characteristics. 


In order to evaluate the weight function in the neighborhood of the 
recollection time origin an analytical integration is obtained when the Bessel 
function is approximated with the inverse of the square root of 2 m times the 
argument, which is equal to 2 rr, /(4 a@ t,). However, this approximation 
imposes a major constraint due to the fact that the argument of the Bessel’s 
function must be greater than 100. For small values of t, the Bessel function 


and its approximation strongly diverge as shown in Fig. 7. 


Since the grid size and the thermal diffusivity are dictated by the 
physical problem, the time step is the only parameter that can be adjusted to 
meet the condition required for the analytical integration of the weight 
function. This means that the time step must be of the order of 0.001 seconds 
which is extremely small compared to the evaporation time. The logical 


alternative is to resort to numerical integration over the initial time step. 


maj 


WEIGHT TEMPORAL CHARACTERISTICS 


yay Gly Awe 
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FIGURE SIX - Weight Temporal Characteristics 
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Lo APPROXIMATION 


CONVERGENT FOR ARGUMENT > 100 


VALUE 


fe) 0.2 0.4 0.6 0.8 1.0 


ARGUMENT = (2 R Ro}/(4@ to) 


FIGURE SEVEN - L, Approximation 
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The discretization of the integrand must have a resolution such that the 
peak of the function can be properly captured. Consequently, the behavior of 
the weight function for various governing parameters versus the recollection 
time becomes essential. In order to describe the weight function for 
combinations of the parameters (Ar,, At,, @), a non dimensional quantity is 
identified as P = Ar,*/(a@ At,). The weight function can now be plotted for 


selected values P within the range of interest, as shown in Fig. 8. 


4.1.3 Matrix Structure Intuitively, the weight function must have its highest 
value when the source coincides with the point of concern and should decay as 
the source moves further away from the point. This behavior is confirmed by 
the numerical values of several test cases. The weights can be arranged in a 
matrix form where each column is composed of all the weights associated with 
the sources affecting one single point. Similarly a row is the collection of 
all the weights associated to one single source affecting the various points 
of the domain. This weight matrix is diagonally dominant and decays on both 
sides of the main diagonal as shown in Table 3. Further, it is important to 
note that, as time elapses the weight matrix values decay rapidly and in fact 

after a few time steps they become negligible. This means that the system 


recollection time extends in the past for a limited number of time steps. 


4.2 Weight Function Validation 


After the preliminary characterization, the weight function must be 


tested in the full range in matrix form to establish its accuracy. 
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WEIGHT TEMPORAL BEHAVIOR 
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WEIGHT 


FIGURE EIGHT - Weight Temporal Behavior 
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TABLE 3 - Weight Matrix Structure 
(Ar, = 100 wp, At, = 1 sec, @ = 2500 p*/sec) 
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TABLE 3 - Weight Matrix Structure, continued 
(Ar, = 100 wp, At, = 1 sec, a@ = 2500 p?/sec) 
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4.2.1 Point source test In order to test the analytical derivation of the 
weights a point source problem is solved and the solution obtained by the 
boundary element method is compared with the available exact solution. The 
exact solution form for a point source is essentially that of the Green's 
function for a point source equation with minor modifications. The 
dimensional units of the exact solution must be adjusted to match those of the 
weight function, namely the dimension of ‘'u’ must be in terms of temperature. 
This is accomplished by first obtaining the solution for the non-dimensional 


governing energy equation: 


—-V'T (20) 


where the temperature is non-dimensionalized with respect to some arbitrary 
temperature difference AT, the length scale is L = (k AT/q"’) then x* = x/L 
and V" = VL. The time is non-dimensionalized with respect L?/a. The point 


source solution for this equation is 


* *. 2 a 
(x - X, ) ERGY = xs 


Tee ( Agee te a (21) 


Finally, by introducing the appropriate dimensional grouping, one obtains 
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2 2 2 
(XE-8X a) te (Yao aet aca ee) 
T=T +S (4a e\eaicaes 4at (25 


The terms ‘S’ is the source strength which is primarily a function of heat 


flux per unit volume and is found to be equal to 


S = AT (qt)? (23) 


By changing to cylindrical coordinates, by setting the source at r, = 0 and by 


imposing that T = 0, the solution on the plane z = 0 is given by 


WE WR Ago 0 ape ae (24) 


Finally, the corresponding flux equation is found by simply differentiating 


the previous equation 


few weg tirareyete? 


ae e 4Gat e 4at (25) 


Zen 


The result obtained imposing this forcing function should now match the one 


obtained from Eq. (24). This is indeed the case as illustrate in Table 4. 
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TO=0.5000 
0.94 0.72 


TO=1.5000 
139 1.07 


TO=2.5000 
Loo 1.39 


TO=3 .5000 
Zale 1.70 


TO=4.5000 
Pal J 2.01 


TO=5.5000 
2.99 2.06 


TO=6.5000 


Seow 2.67 


TO=7.5000 
3278 SO. 


TO=8 .5000 
4.08 Jee 


TO=9.5000 
4.08 3-26 


EXACT 
3.9787 
3.2575 
2.1836 
1.1984 
0.5385 
0.1981 
0.0597 
0.0147 


TABLE 4 - Weight Function Validation (Ar, = 100 w, At, = 1 sec, 
S=1.9 x 108 » °C @ z, = - 200 pz) 
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4.2.2 Comparison with aluminum results 


The next logical step would be to use the weight function to obtain an 
analytical solution for the aluminum experiment and check the results 

against the existing data base. The aluminum problem is a desirable one 

since the solution for the solid thermal behavior is uncoupled from the 
droplet evaporation process. This is due to the fact that for materials with 
thermal conductivity the temperature at the solid-liquid interface remains 
uniform and constant during the evaporative transient. Therefore no 
information is needed on the liquid side and the solid side weight function 


will suffice for the solution. 


The semi-infinite surface of the aluminum can be divided into two 
regions: a) wetted region (under the droplet), where the temperature is known 
and the b) exposed region, where the flux is known. At this point an 
appropriate transformation of the temperature must be chosen in order to 
eliminate the volume integral at the initial time in Eq. (8) as illustrated in 
Eq. (9). Consider the following expression of the temperature as a function 


of depth 


Ge ae (26) 


and define u as 


u=T- T, (27) 


oy 


then the flux at the exposed solid surface is 


du 
az 


--2u (28) 


note that T, is the initial temperature at each point of the solid domain. 
Equation (28) is written at the surface of the domain where the boundary 


condition must be imposed. 


To illustrate the solution scheme for the semi-infinite solid, a domain 
on the solid surface is chosen. The domain consists of 5 points of which 3 are 
under the droplet and 2 are on the exposed surface. The equation written for 


these five points are 


u, = Wi, f, + wf, +w,,f; + w,,f, + wf + r, 
u, = wo,f, + wo,f, + w23f, + w,,f, + w,5f, + rz 
Ug Wey Lee WT, Pa et wy sf, ew, SEY wet re (29) 
uy * Wot, “tow fst Ws f, Pty, f, + west. + ir, 


u, = ws,f, + Wof, + ws3f, > Wo, t, + woof. tare 


The unknowns are the fluxes applied on the wetted region and the temperatures 
on the exposed surface, namely: i eee od RES u,. Regrouping the unknowns 


on the left hand side, one obtains 
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Wi,f, + W,2f, + w,3f, mm We Xin 


Woif, + w.,f, + w23f, “=U, - T2 - W,,f, - Wos5fs 
W3if, + w3,f, + w,f, —u, - 13 - W3,f, - wys5fs (30) 
wif, + wj2f + wf; -u- Ty Coa ays eyes 
Ws5if, + Ws,f, + we,f,; - us = - fy - Ws5,f, - Wests 


Wii Wy2 W3 0 0 fy u - T, - Wf, - wif 
Wei We2 W23 9 O f, uz - Tz - w,f, - wf 
W31 W532 W33 0 O fz | = | 43 - 53 - W3,f, - W35f (31) 
We. Wer Wz “lL O u, Se eae Pee 
: Wsi Ws2 Ws53 9 -1 Us - Ts - W5,f, - wesfs 


The solution to this matrix expression is obtained by inversion and typical 
results are shown in Table 5. These results are in good agreement with the 
ones obtained with the finite difference method described by diMarzo 1986a. 
This first solution is used to validate the boundary element technique which 
will be extended to the more complex coupled solution for the low thermal 


conductivity solid case. 


4.3 Codes 
The codes used in this portion of the study are enclosed for reference. 
In particular the weight matrix generator, the point source validation and the 


semi-infinite solid solution used for the aluminum case are included. 


T= 0.50 
SOLUTION X 
-1.445 -1.662 -4.124 
SOLUTION U OR T-TO 
-6.000 -6.000 -6.000 


T= 1.50 
SOLUTION X 
-1.368 -1.575 -3.924 -0.828 
SOLUTION U OR T-TO 
-6.000 -6.000 -6.000 -0.828 


T= 2.50 
SOLUTION X 
-1.336 -1.539 =-3.838 -2.404 
SOLUTION U OR T-TO 
-6.000 -6.000 -6.000 -2.404 


T= 3.50 
SOLUTION X 
~1. 5195-12519 esas ~2.449 -0.983 
SOLUTION U OR T-TO 
-6.000 -6.000 -6.000 -2.449 -0:983 


T= 4.50 
SOLUTION X 
-1.307 “=1.505 =S2756 “2.479 
SOLUTION U OR T-TO 
-6.000 -6.000 -6.000 2.479 


T= 5.00 
SOLUTION X 
~1 29801, A96--3 9152 -2.501 -1.052 
SOLUTION U OR T-TO 
-6.000 -6.000 -6.000 2.501 -1.0352 


T=" 6.50 
SOLUTION X 
“1.292 ~-1.488 1=3.714 =2°518 -1.075 
SOLUTION U OR T-TO 
-6.000 -6.000 -6.000 27916 =1.079 


T= 7.50 
SOLUTION X 
-1.287 -1.482 -3.699 -1.094 
SOLUTION U OR T-TO 
-6.000 -6.000 -6.000 1-30.94 


TABLE 4 - Solution of the Aluminum Case; U is the temperature drop, 
T is the time (Ar, = 1 mm, At, = 1 sec, a = 45 mm*/sec) 


4.3.1 Weight matrix generator 
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PROGRAM WTI 

REAL W(10,10) 

COMMON /CONST/ ALFA, PI 

COMMON /GRID/ DRO,DT0O,NP 

COMMON /INTG/ ATO,BTO,AERR, RERR 


PI=3.1415 

ALFA=2500.0 

READ (3, *) DRO,DT0,NP, TN 
READ (3, *) ATO, BTO, AERR, RERR 
TS=TN/DTO0+0.01 

WRITE (9,21) DRO,DTO,NP 


DO 80 TK=1.0,TS 
TO=(TK-0.5) *DTO 
WRITE (9, 61) TO 


CALL WEIGHT (TO, W) 


WRITE (9, 81) ((W(I,J),J=1,NP) , I=1,NP) 


CONTINUE 


FORMAT (2X, ’DRO=’,F7.2,2X,’DT0=’,F5.2,2X, 'NP=’, 14) 


FORMAT (/,5X,‘’TO=’,F6. 4) 
FORMAT (8 (2X, F6.2) ) 

STOP 

END 


SUBROUTINE WEIGHT (TO, W) 


REAL W(10,10) 
REAL DCADRE, MMBSI0,ERF 
EXTERNAL F 


COMMON /CONST/ ALFA, PI 

COMMON /GRID/ DRO,DTO,NP 

COMMON /INTG/ ATO,BTO,AERR, RERR 
COMMON /RRO/ R,RO 


ARG1=DR0/SQRT (16*ALFA*TO) 
DO 50 I=1,NP 

DO 50 J=1,NP 

R= (FLOAT (I) -0.5) *DRO 

RO= (FLOAT (J) -0.5) *DRO 
ARG3= (R-RO) **2/ (4*ALFA*TO) 
ARG4=0 .5*R*RO/ (ALFA*TO) 


SET UP (W) FOR TO < DTO 


IF (TO .LT.DTO) THEN 


W(I,J) =DCADRE (F, ATO, BTO, AERR, RERR, ERROR, IER) 


END IF 


SET UP (W) FOR TO > DTO 


50 


IF (T0.GE.DT0) THEN 
IF (I.NE.J) THEN 


W(I,J)=(1.0/SQRT(4.0*PI*ALFA) ) *RO* (TO**(-1.5)) * 
MMBSTI0 (2, ARG4, IER) *EXP (-ARG3) *DRO*DTO 


ELSE 


W(I,J)=(R/TO) *MMBSIO (2, ARG4, IER) *ERF (ARG1) *DTO 


END IF 
END IF 


CONTINUE 
RETURN 
END 


FUNCTION F (TO) 
REAL MMBSI0O, ERF 


WTI00010 
WTI00020 
WTI00030 
WTI00040 
WTI00050 
WTI00060 
WT1I00070 
WTI00080 
WTI00090 
WTI00100 
WTI00110 
WTI00120 
WTI00130 
WTI00140 
WTI00150 
WTI00160 
WTI00170 
WTI00180 
WTI00190 
WTI00200 
WTI00210 
WTI00220 
WTI00230 
WTI00240 
WTI00250 
WTI00260 
WT1I00270 
WTI00280 
WTI00290 
WTI00300 
WTI00310 
WTI00320 
WTI00330 
WTI00340 
WTI00350 
WTI00360 
WTI00370 
WTI00380 
WTI00390 
WTI00400 
WTI00410 
WTI00420 
WTI00430 
WTI00440 
WTI00450 
WTI00460 
WTI00470 
WTI00480 
WTI00490 
WTI00500 
WTI00510 
WTI00520 
WTI00530 
WTI00540 
WTI00550 
WTI00560 
WTI00570 
WTI00580 
WTI00590 
WTI00600 
WTI00610 
WTI00620 
WTI00630 
WTI00640 
WT 

WI 42 


COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NP 
COMMON /RRO/ R,RO 


ARG1=DR0/SQRT (16*ALFA*TO) 

ARG3=(R-RO) **2/ (4*ALFA*TO) 

ARG4=0 .5*R*RO/ (ALFA*TO) 

E=0.0 

IF (ARG3.LE.150.0) E=EXP (-ARG3) 

IF (R.NE.RO) THEN 

F=(1.0/SQRT(4.0*PI*ALFA) ) *RO* (TO**(-1.5))* 
MMBSI0 (2, ARG4, IER) *E*DRO 

ELSE 

F=(R/TO) *MMBSIO (2, ARG4, IER) *ERF (ARG1) 

END IF 

RETURN 

END 


WTI00670 
WTI00680 
WTI00690 
WTI00700 
WTI00710 
WTI00720 
WTI00730 
WTI00740 
WTI00750 
WTI00760 
WTI00770 
WTI00780 
WTI00790 
WTI00800 
WTI00810 
WTI00820 
WTI00830 
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PROGRAM WT 

REAL W(10,10) 

COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NP 


PI=3.1415 

ALFA=2500.0 

READ (3, *)DRO,DTO,NP, TN 
TS=TN/DTO+0.01 

WRITE (8,21)DRO,DTO,NP 


DO 80 TK=1.0,TS 
TO=(TK-0.5) *DTO 
WRITE (8, 61) TO 


CALL WEIGHT (TO, W) 
WRITE (8,81) ((W(I,J) ,J=1,NP) , I=1,NP) 
CONTINUE 


FORMAT (2X, ’DRO=’ ,F7.2,2X, 'DT0=",F5.2,2X, ’NP=’ ,T4) 
FORMAT (/,5X,‘’TO=’,F6. 4) 

FORMAT (8 (2X, F 6.3) ) 

STOP 

END 


SUBROUTINE WEIGHT (TO, W) 


REAL W(10,10) 

REAL MMDEI,ERF,MMBSIO 
COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NP 


9 SET UP (W) FOR TO < DTO 


20 


IF (TO.LT.DTO) THEN 
ARG1=DR0/SQRT (16*ALFA*DTO) 
ARG2=ARG1**2 


DO 20 I=1,NP 

DO 20 J=1,NP 

R= (FLOAT {I) -0.5) *DRO 

RO= (FLOAT (J) -0.5) *DRO 

ARG4= (R-RO) **2/ (4*ALFA*DTO) 

IF (I.NE.J) THEN : 

W(I,J)=1.0/ (2.0*PI) *SQRT(RO/R) *MMDEI (2, ARG4, IER) *DRO 

ELSE 

W(I,J)=2.0*SQRT (ALFA/PI) *SQRT (DTO) *ERF (ARG1) + 
DRO/ (2.0*PI) *MMDEI (2, ARG2, IER) 

END IF 

CONTINUE 

END IF 


te SET UP (W) FOR TO > DTO 


IF (TO.GE.DTO) THEN 

ARG1=DR0/SQRT (16*ALFA*TO) 

DO 50 I=1i,NP 

DO 50 J=1,NP 

R= (FLOAT (I) -0.5) *DRO 

RO= (FLOAT (J) -0.5) *DRO 

ARG3= (R-RO) **2/ (4*ALFA*TO) 

ARG4=0 .5*R*RO/ (ALFA*TO) 

IF (I.NE.J) THEN 

W(I,J)=(1.0/SQRT(4.0*PI*ALFA) ) *RO* (TO**(-1.5))* 
MMBSI0 (2, ARG4, IER) *EXP (-ARG3) *DRO*DTO 

ELSE 

W(I,J)=(R/TO) *MMBSIO (2, ARG4, IER) *ERF (ARG1) *DTO 

END IF 


WEI00010 
WEI00020 
WEI0N0030 
WEI00040 
WEI00050 
WEI00060 
WEI00070 
WEI00080 
WEI00090 
WEI00100 
WEI00110 
WEI00120 
WEI00130 
WEIO0140 
WEI00150 
WEI00160 
WEI00170 
WEI00180 
WEI00190 
WEI00200 
WEI00210 
WEI00220 
WEI00230 
WEI0N0240 
WEI00250 
WEI00260 
WEI00270 
WEI0N0280 
WEI00290 
WEI00300 
WEI00310 
WEI00320 
WEI00330 
WEI00340 
WEI00350 
WEI00360 
WEI00370 
WEI0N0380 
WEI00390 
WEI00400 
WEI00410 
WEI00420 
WEI00430 
WEI00440 
WEI00450 
WEI0N0460 
WEI00470 
WEI00480 
WEI00490 
WEI00500 
WEI00510 
WEI0N0520 
WEI00530 
WEI00540 
WEI00550 
WEI00560 
WEI00570 
WEIO00S580 
WEIO0590 
WEI00600 
WEI00610 
WEI00620 
WEI00630 
WETIANEAN 
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CONTINUE 
END IF 
RETURN 
END 


WEI00670 
WEI00680 
WEI00690 
WEI00700 


45 


‘> 
» 


f 


v Fi 
" DAU 2S, 2.53 


1 ¥ 


we. 


174s on 
Le | +4 

s 
Bidet 
Wide ioe 
ate 9% 


aS 


* 


PROGRAM MAIN 


POINT SOURCE TEST FOR THE WEIGHT FUNCTION 
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REAL W(10,10),F (10) ,WF(10),R(10) ,U(10) 
COMMON /CONST/ ALFA, PI 

COMMON /GRID/ DRO,DTO,NP,TN 

COMMON /INTG/ ATO,BTO,AERR, RERR 
COMMON /SOURCE/ S,Z0 


PI=3.1415 

ALFA=2 .5E+03 

READ (3, *)DRO,DT0,NP,TN 
READ (3, *) ATO, BTO, AERR, RERR 
READ (3, *) RN, Z0 

WRITE (10,21)DRO,DTO,NP,TN 
WRITE (10, 31) RN, ZO 


TS=TN/DT0+0.01 
S=(RN/SQRT (12.0) ) **3 


DO 20 I=1,NP 
R(I)=0.0 

WF (I) #0 .0 
CONTINUE 


DO 80 TK=1.0,TS 


TO= (TK=0'.15):40TC 
WRITE (10, 41) TO 


SET UP THE WEIGHT MATRIX 


CALL WEIGHT (TO, W) 
WRITE (10,61) ((W(I,J),J=1,NP) ,I=1,NP) 


COMPUTE THE FLUXES 


CALL FLUX(TO,F) 
WRITE (10, 61) (F (J), J=1,NP) 


“SUMMATION IN SPACE 


CALL MULT (NP,W,F, WF) 
WRITE (10, 61) (WF (J), J=1,NP) 


SUMMATION IN TIME 


80 


CALL ADD (NP, WF, R) 
WRITE (10,61) (R(J) , J=1,NP) 


CONTINUE 


COMPUTE THE EXACT SOLUTION 


WRITE (10, *) 

WRITE (10,*) 7 EXACT 
CALL EXACT (U) 

WRITE (10,91) (U(J),R(J) , J=1,NP) 


NUMERICAL’ 


FORMAT (2X; DRO=*,F7.2,24,, D100" ,F5.2,/,2%,.. NP=".,14, 


bag. TN Ege) 
FORMAT (2X, 'RN#=’, F7.2,2X,'Z0=#’ ,F7.2) 
FORMAT (/,5X,’TO=’,F6.4) 
FORMAT (8 (2X, F6.2) ) 
FORMAT (11X,F7.4,14X,F7. 4) 
STOP 
END 


SUBROUTINE WEIGHT (TO,W) 


REAL W(10,10) 
REAL DCADRE,MMBSI0,ERF 
EXTERNAL F 


COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NP,TN 
COMMON /INTG/ ATO,BTO,AERR, RERR 


TINO0010 
TINO0020 
TINO0030 
TINO0040 
TINO0050 
TINO0060 
TINO0070 
TINO0080 
TINO0090 
TINO0100 
TINO0110 
TINO0120 
TINO0130 
TINO0140 
TINO0150 
TINO0160 
TINO0170 
TINO0180 
TINO0190 
TINO0200 
TINO0210 
TINOQ0220 
TIN00230 
TINO0240 
TINO0250 
TINO0260 
TINO0270 
TINO0280 
TINO0290 
TINO0300 
TINO0310 
TINO0320 
TINO00330 
TINO0340. 
TINO0350 
TINO0360 
TINO0370 
TINO0380 
TIN00390 
TINO0400 
TINO0410 
TINO0420 
TINO0430 
TINO0440 
TINO0450 
TINO0460 
TINO0470 
TINO0480 
TINO0490 
TINOOSO0 
TINO0510 
TINOO0S520 
TINO0530 
TINO0540 
TINOO0550 
TINOOS560 
TINO0570 
TINO0S580 
TINO0590 
TINO0600 
TINO0610 
TINO0620 
TINO0630 
TINON640 
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COMMON /RRO/ R,RO 


ARG1=DR0/SORT (16*ALFA*TO) 
DO 50 I=1,NP 

DO 50 J=1,NP 

R= (FLOAT (I) -0.5) *DRO 

RO= (FLOAT (J) -9.5) *DRO 
ARG3=(R-RO) **2/ (4*ALFA*TO) 
ARG4=0 .5*R*RO/ (ALFA*TO) 


SET UP (W) FOR TO < DTO 


IF (TO.LT.DTO) THEN 
W(I,J) =DCADRE (F, ATO, BTO, AERR, RERR, ERROR, IER) 
END IF 


SET UP (W) FOR TO > DTO 


50 


50 


IF (TO .GE.DT0O) THEN 

IF (I.NE.J) THEN 

W(I,J)#(1.0/SQRT (4.0*PI*ALFA) ) *RO* (TO** (-1.5)) * 
MMBSI0 (2, ARG4, IER) *EXP (~ARG3) *DRO*DTO 

ELSE 

W(I,J)=(R/TO) *MMBSIO (2, ARG4, IER) *ERF (ARG1) *DTO 

END IF 

END IF 


CONTINUE 
RETURN 
END 


FUNCTION F (TO) 

REAL MMBSI0,ERF 

COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NP,TN 
COMMON /RRO/ R,RO 


ARG1=DR0/SQRT (16*ALFA*TO) 
ARG3=(R-RO) **2/ (4*ALFA*TO) 
ARG4=0 .5*R*RO/ (ALFA*TO) 

E=0.0 

IF (ARG3.LE.150.0) E=EXP (-ARG3) 


IF (R.NE.RQ) THEN 

F=(1.0/SQRT(4.0*PI*ALFA) ) *RO* (TO**(-1.5)) * 
MMBSI0 (2, ARG4, IER) *E*DRO 

ELSE 

F=(R/TO) *MMBSIO(2,ARG4, IER) *ERF (ARG1) 

END IF 


SUBROUTINE FLUX(TO,F) 


REAL F(10) 

COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NP,TN 
COMMON /SOURCE/ S,2Z0 


T=TN-TO 

DO 50 J=1,NP 

RO= (FLOAT (J) -0.5) *DRO 

C=4.0*ALFA*T 

F(J)=5*(C** (-1.5)') * (2*20/C) EXE = (20**2) /C) See e (Roe eer) 
CONTINUE 

RETURN 

END 


SUBROUTINE MULT (NP,W,F, WF) 


TINO0670 
TINO0680 
TINO0690 
TINOO700 
TINOO710 
TINOO720 
TINOO730 
TINOO740 
TINO075S0 
TINOO760 
TINO0770 
TINOO780 
TINO0790 
TINOO800 
TINO0810 
TINO0820 
TINO0830 
TINO0840 
TINOO8S0 
TINOO860 
TINO0870 
TINO0880 
TINO0890 
TINOO900 
TINO0910 
TINOO920 
TINO0930 
TINO0940 
TINOO950 
TINOO960 
TINO0970 
TINOO980 
TINO0990 
TINO1000 
TINO1010 
TINO1020 
TINO1030 
TINO1040 
TINO1050 
TINO1060 
TINO1070 
TINO1080 
TINO1090 
TINO1100 
TINO1110 
TINO1120 
TINO1130 
TINO1140 
TINO1150 
TINO1160 
TINO1170 
TINO1180 
TINO1190 
TINO1200 
TINO1210 
TINO1220 
TINO1230 
TINO1240 
TINO1250 
TINO1260 
TINO1270 
TINO1280 
TINO1290 
TINO1300 
Tis 


50 


50 


50 


REAL W(10,10),F (10) ,WF (10) 


DO 50 I=1,NP 

WF (I) =0.0 

DO 50 K=1,NP 

WE (I) =WF (I) +W(I,K) *F (K) 
CONTINUE 

RETURN 

END 


SUBROUTINE ADD (NP, WF, R) 
REAL WF(10),R(10) 


DO 50 I=1,NP 
R(I) #R (I) +WE (T) 
CONTINUE 
RETURN 

END 


SUBROUTINE EXACT (U) 

REAL U(10) 

COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NP,TN 
COMMON /SOURCE/ S,2Z0 


DO 50 I=1,NP 
R= (FLOAT (I) -0.5) *DRO 
C=4.0*ALFA*TN 


Utl mor (ee (=) Exe (=(20%*2) 7 C) XEXP (- (R*¥*2):/C) 


CONTINUE 
RETURN 
END 


TINO1330 
TINO1340 
TINO1350 
TINO1360 
TINO1370 
TINO1380 
TINO1390 
TINO0O1400 
TINO1410 
TINO1420 
TINO1430 
TINO1440 
TINO1450 
TINO1460 
TINO1470 
TINO1480 
TINO1490 
TINO1500 
TINO1510 
TINO1520 
TINO1530 
TINO1540 
TINO1550 
TINO1560 
TINO1570 
TINO1580 
TINO1590 
TINO1600 
TINO1610 
TINO1620 
TINO1630 
TINOQ1640 
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PROGRAM ALI ALI00010 


* ALI00020 
REAL W(50,50),F(50),FTIM(50, 60) ,WE (50) ,R(50) ,U(50) ALI00030 
REAL WO (50,50) ,AINV(50,50),A(50,50) , WKAREA(50) ALI00040 
REAL X(50),B(50) ALI00050 
COMMON /CONST/ ALFA,PI ALI00060 
COMMON /GRID/ DRO,DT0,NDP,NP,TN ALI00070 
COMMON /INTG/ ATO,BTO,AERR, RERR ALI00080 
* ALI00090 
PI=3.1415 ALI00100 
ALFA=45 .0 ALI00110 
HK=1.1E-04 ALI00120 
READ (2, *)DRO,DT0,NDP,NP, TN ALI00130 
READ (2, *) ATO, BTO, AERR, RERR ALI00140 
WRITE (20,21) DRO,DTO,NDP,NP, TN ALI00150 
* ALI00160 
IEND=INT (TN/DTO) ALI00170 
NT=INT (TN) ALI00180 
* ALI00190 
* INITIALIZE ALI00200 
DO 10 I=1,NP ALI00210 
R(I)=0.0 ALI00220 
WE (I) =0.0 ALI00230 
F(I)=0.0 ALI00240 
U(I)=0.0 ALI00250 
10 CONTINUE AL1I00260 
DO 20 I=1,NDP AL1I00270 
U(I)=-6.0 ALI00280 
20 CONTINUE ALI00290 
* ALI00300 
* COMPUTE THE LARGEST WEIGHT & STORE IT ALI00310 
TO=0.5*DTO ALI00320 
CALL WEIGHT (TO, W) ALI00330 
DO 50 I=1,NP ALI00340 
DO 50 J=1,NP ALI00350 
WO (I,J) =W(I,J) ALI00360 
50 CONTINUE ALI00370 
WRITE (20,81) ((WO(I,J),J=1,NP) , I=1,NP) ALI00380 
x ALI00390 
DO 200 KT=1, IEND ALI00400 
T= (FLOAT (KT) -0.5) *DTO ALI00410 
WRITE (20,61)T ALI00420 
* FIND MOST RECENT FLUX ALI00430 
CALL SETUP (W0O,U,F,R,A,B) ALI00440 
* RESET R(I)=0.0 ALI00450 
DO 60 I=1,NP ALI00460 
R(I)=0.0 ALI00470 
60 CONTINUE ALI00480 
* ALI00490 
IF (KT.EQ.1) THEN ALI00500 
CALL LINV1F(A,NP,50,AINV, IDGT, WKAREA, IER) ALI00510 
ENDIF ALI00520 
* ALI00530 
CALL MULT (NP,NP,AINV,B, X) ALI00540 
* ALI00550 
WRITE (20,*)’ SOLUTION X’ ALI00560 
WRITE (20, 91) (X(I) , I=1, NP) ALI00570 
* UPDATE THE SOLUTION ALI00580 
DO 70 I=1,NP ALI00590 
IF (I. LE.NDP) F (I) =X (I) ALI00600 
IF (I.GT.NDP) THEN ALI00610 
U(I) =X (I) ALI00620 
F (I) =-HK*U (I) ALI00630 
ENDIF Bit Gee ** 

70 CONTINUE ALI 


* STORE FLUXES FOR EACH TIME T ALI 51 


* 


80 


DO 80 I=1,NP 
FTIM(I,KT) =F (I) 
CONTINUE 


SET UP R FOR NEXT CALCULATION 


DO 150 IT=1,KT 


SELECT PROPER WEIGHTS 


TO=FLOAT (KT-IT)+1.5 
CALL WEIGHT (TO, W) 


PICK UP CORRESPONDING FLUX 


100 


DO 100 I=1,NP 
F(I)=FTIM(I, IT) 
CONTINUE 


SUMMATION IN SPACE 


CALL MULT (NP,NP,W,F, WE) 


SUMMATION IN TIME 


150 


60 


80 


100 


CALL ADD (NP, WF, R) 


CONTINUE 
WRITE (20,*)’ SOLUTION U OR T-TO’ 
WRITE (20,91) (U(I) , I=1,NP) 


CONTINUE 

WRITE (20, *) ‘FLUXES’ 

WRITE (20,81) ((FTIM(I, KT), I=1,NP) ,KT=1,NT) 

FORMAT (2X, ’DRO=’,F7.2,2X,’DT0=’,F5.2,2X,/,’NDP=’ , 14, 2X, 
'DP=’,14,2X,’TN=’,F5.2) 

FORMAT (/,5X,’T=',F6.2) 

FORMAT (9 (1X, F6.3) ) 

FORMAT (9 (1X, F6.3) ) 

STOP 

END 


SUBROUTINE SETUP (W0,U,F,R,A,B) 

REAL WO (50,50),F(50) ,WE(50),R(50) ,U(50) 
REAL A(50,50),WE(50, 45) ,FE(45),¥(50),B(50) 
COMMON /GRID/ DRO,DTO,NDP,NP,TN 


DO 60 I=1,NP 

DO 60 J=1,NP 

A(I,J)=0.0 
IF(J.LE.NDP)A(I,J) =W0 (I,J) 
IF(J.GT.NDP .AND. I.EQ.J)A(I,J)=-1 
CONTINUE 

NEP=NP-NDP 

DO 80 I=1,NP 

DO 80 J=1,NEP 

WE (I, J) =WO (I, J+NDP) 

FE (J) =F (J+NDP) 

CONTINUE 

CALL MULT (NP, NEP, WE, FE, Y) 
C=1.0 

DO 100 I=1,NP 

IF (I.GT.NDP) C=0.0 

Bi D)eC*U (1) =R (1) —2 (1) 
CONTINUE 

RETURN 

END 


SUBROUTINE MULT (M,N,W,F, WF) 
REAL W(50,50),F (50) ,WE (50) 


DO 50 I=1,M 
WF (I) =0.0 


ALI00670 
ALI00680 
ALI00690 
ALI00700 
ALI00710 
ALI00720 
ALI00730 
ALI00740 
ALI00750 
ALI00760 
ALI00770 
ALI00780 
ALI00790 
ALI00800 
ALI00810 
ALI00820 
ALI00830 
ALI00840 
ALI00850 
ALI00860 
ALI00870 
ALI00880 
ALI00890 
ALI00900 
ALI00910 
ALI00920 
ALI00930 
ALI00940 
ALI00950 
ALI00960 
ALI00970 
ALI00980 
ALI00990 
ALI01000 
ALI01010 
ALI01020 
ALI01030 
ALI01040 
ALI01050 
ALI01060 
ALI01070 
ALI01080 
ALI01090 
ALI01100 
ALI01110 
ALI01120 
ALI01130 
ALI01140 
ALI01150 
ALI01160 
ALI01170 
ALI01180 
ALI01190 
ALI01200 
ALI01210 
ALI01220 
ALI01230 
ALI01240 
ALI01250 
ALI01260 
ALI01270 
ALI01280 
ALI01290 
ALIN1320NN 


AL 
ys a haben 


50 


50 


DO 50 K=1,N 

WF (I) =WF (I) +W(I,K) *F (K) 
CONTINUE 

RETURN 

END 


SUBROUTINE ADD (NP, WF, R) 
REAL WE (50),R(50) 


DO 50 I=1,NP 
R(I) =R(1) +WF (TI) 
CONTINUE 
RETURN 

END 


SUBROUTINE WEIGHT (TO, W) 


REAL W(50,50) 
REAL DCADRE, MMBSI0,ERF 
EXTERNAL F 


COMMON /CONST/ ALFA, PI 

COMMON /GRID/ DRO,DTO,NDP,NP,TN 
COMMON /INTG/ ATO,BT0O,AERR, RERR 
COMMON /RRO/ R,RO 


ARG1=DR0/SQRT (16*ALFA*TO) 
DO 50 I=1,NP 

DO 50 J=1,NP 

R= (FLOAT (I) -0.5) *DRO 

RO= (FLOAT (J) -0.5) *DRO 
ARG3= (R-RO) **2/ (4*ALFA*TO) 
ARG4=0 .5*R*RO/ (ALFA*TO) 


SET UP (W) FOR TO < DTO 


-IF(TO.LT.DTO) THEN 
W(1I,J) =DCADRE (F, ATO, BTO, AERR, RERR, ERROR, IER) 
END IF 


SET UP (W) FOR TO > DTO 


50 


IF (TO.GE.DTO) THEN 

IF (I.NE.J) THEN 

W(I,J)=(1.0/SQRT (4.0*PI*ALFA) ) *RO* (TO** (-1.5)) * 
MMBSI0 (2, ARG4, IER) *EXP (-ARG3) *DRO*DTO 

ELSE 

W(I,J)=(R/TO) *MMBSTIO (2, ARG4, IER) *ERF (ARG1) *DTO 

END IF 

END -IF 


CONTINUE 
RETURN 
END 


FUNCTION F(TO) 

REAL MMBSIO,ERF 

“COMMON /CONST/ ALFA, PI 

COMMON /GRID/ DRO,DT0,NDP,NP,TN 
COMMON /RRO/ R,RO 


ARG1=DR0/SQRT (16*ALFA*TO) 
ARG3=(R-RO) **2/ (4*ALFA*TO) 
ARG4=0 .5*R*RO/ (ALFA*TO) 

E=0 .0 

IF (ARG3.LE.150.0) E=EXP (-ARG3) 


IF (R.NE.RO) THEN 
F=(1.0/SQRT(4.0*PI*ALFA) ) *RO* (TO** (-1.5) ) * 


ALI01330 
ALI01340 
ALI01350 
ALI01360 
ALI01370 
ALI01380 
ALI01390 
ALI01400 
ALI01410 
ALI01420 
ALI01430 
ALI01440 
ALI01450 
ALI01460 
ALI01470 
ALI01480 
ALI01490 
ALI01500 
ALI01510 
ALI01520 
ALI01530 
ALI01540 
ALI01550 
ALI01560 
ALI01570 
ALI01580 
ALI01590 
ALI01600 
ALI01610 
ALI01620 
ALI01630 
ALI01640 
ALI01650 
ALI01660 
ALI01670 
ALI01680 
ALI01690 
ALI01700 
ALI01710 
ALIQ1720 
ALI01730 
ALI01740 
ALI01750 
ALI01760 
ALI01770 
ALI01780 
ALI01790 
ALI01800 
ALI01810 
ALI01820 
ALI01830 
ALI01840 
ALI01850 
ALI01860 
ALI01870 
ALI01880 
ALI01890 
ALI01900 
ALI01910 
ALI01920 
ALI01930 
ALI01940 
ALI01950 
ALIQ1Q4N 
AL: 

AL. 53 


MMBSI0 (2, ARG4, IER) *E*DRO 
ELSE 
F=(R/TO) *MMBSIO (2, ARG4, IER) *ERF (ARG1) 
END IF 
RETURN 
END 


ALI01990 
ALI02000 
ALI02010 
ALI02020 
ALI02030 
ALI02040 
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a CONCLUSION AND FUTURE WORK 

The preliminary results obtained with the boundary element method are 
very promising. This analytical task will be carried on in the next grant 
period. The coupled solution for the low thermal conductivity solid case will 
be obtained and compared with experimental measurements of the solid surface 
temperature via infrared thermography. A new set of experiments on a quartz 
semi-infinite solid are under way and initial efforts to set-up the infrared 
thermography apparatus have been made. Another area of research that will be 
developed in the next grant period is involved with the generation of small 
droplets (order of one microliter). These droplet are in the same size range 


of those generated by commercial sprinklers. 


The solids analyzed in this initial portion of the research program are 
heated from below by an electrical hot plate. In the actual fire application 
the heat is transmitted to the water droplet and to the surface by radiation 
from above. A new experimental set-up where both the droplet and the surface 
will receive radiative heat from gas fired panels is planned for the end of 
the next grant period. A solid block with embedded thermocouples will also be 
tested in this radiant heat set-up to validate the modelling of the solid 


thermal behavior for internal points. 


=P) 
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